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The impetus for this study came from the need to accurately predict the 
performance of propellers on V/STOL aircraft operating in the static 
condition. Small errors in thrust estimation are easily magnified into 
large errors in payload estimation. At the start of this study, it was 
felt that classical propeller analyses as well as some more recent numerical 
analyses methods did not adequately predict the static performance. 

The classical vortex theory analyses for propellers aie based on the 
physical situation of having the propeller advance at stn'ne finite forward 
velocity. In this theory each blade is modelled as a straight bound vortex 
filament and the wake behind each blade is represented by a force-free 
vortex sheet. For a lightly loaded optimum propeller, Betz 1 showed that 
the geometry of each trailing vortex sheet is that of an uncontracted 
helical surface which is aligned with the resultant velocity in the slip- 
stream. Goldstein ' 1 was able to calculate the performance for the lightly 
loaded optimum propeller and Theodor sen 3 later extended Goldstein's analysis 
to predict the performance of moderately "loaded propellers, still retaining- 
the true helical surfaces as the model of the trailing vortex sheets. When 
the classical methods are applied to the statically thrusting propeller, 
the predictions tend to be optimistic. 

In the actual static propeller case, experimental evidence shows that 
the wake has a high degree of contraction, the vortex sheets near the tip 
tend to roll up into strong discreet vortices, and the inner part of each 
sheet tends to be distorted from the classical helical model. 

In an attempt to more accurately model the static wake two other 
distinct approaches have generally been tried — the free-wake methods and 
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the prescribed-wake methods. The free-vortex approach is exemplified by 
Erickson and Ordway 1 *. Their work if based on vortex theory in which they 
fix the wake contraction by means of heavily loaded actuator disk theory. 

The force-free condition for the trailing vortex sheets is then obtained by 
iterating on the pitch of these sheets. Characteristically this approach 
requires a large number of iterations and the final results ate somewhat 
dependent upon the original chosen form of the vortex sheets. 

The prescribed-wake analysis, is a semi-empirical approach exemplified 
by the works of Landgrebe 5 and Ladden 6 . In this approach the wakes 
observed empirically are modelled and used to determine the induced velocity 
picture at the propeller blades. This kind of approach has the distinct 
advantage of using little computer time. These methods are, however, 
somewhat dependent upon the availability of experimental wake data. 

—In view of the need to eliminate. any assumptions or empirical restric- 
tions regarding wake shape, the main thrust of the present Investigation 
has been to generate a wake without these restrictions and, thus, be able 
to calculate the induced flow at the propeller blades. To do this it is 
noted that the inflow is known exactly at one Instant of time for any 
propeller; namely, at the instant of start of the propeller motion. Since 
no wake exists at this instant, the inflow is entirely determined by the 
blade motion. As the motion progresses, the wake is deposited and deforms 
continuously under its own self-induced effects until a final shape such 
as observed in Reference 5 is established. This means that the inflow and " 
therefore the loading change continuously until the final wake is established 
and a steady state performance is reached. Essentially, the wake formation 
is treated as an initial condition problem in time. Such a formulation 


4 


implies an unsteady aerodynamic analysis for propellers similar to the 
Wagner problem of fixed-wing aerodynamics. 

This report primarily summarizes the efforts toward this end. The 
principal recording of this work is in Reference 7. Prior to embarking 
on this study it was felt appropriate to develop a more standard 
performance calculation procedure to be used as a point of reference. The 
result was the somewhat modified prescribed wake procedure reported in 
Reference 8. Finally, in an effort to clear up some details of the program 
in Reference 7, particularly in regard to distortions of the vortex filaments 
that occur in the wake and the specification of core sizes for these 
filaments, Reference 9 was written. 

A Reference Static Performance Method 

In Reference 8, Miller reported on the development of a simple numerical 
method to rapidly predict the. static performance of propellers. The wake 
model used in this development is essentially that of Gray 10 who quantified 
the geometry of the tip vortex as a function of the performance of the 
propeller. However, Miller represented the rolled up tip vortices as a 
series of vortex rings whose position was consistent with the quantification 
proposed by Gray. In effect then the rings are used to calculate the axial 
component of induced velocity that would be produced by the rolled up tip 
vortices. The tangential component of induced velocity is derived by an 
analysis similar to that used in the classical vortex theory where the 
entire vortex sheet that is originally shed from each blade is taken into 
account. Finallv, normality is imposed with regard to finding the additional 
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axial induced velocity that corresponds with the tangential component just 
calculated. This means that the vortex ring axial induced component is 
divorced from the normality consideration. 

In order to evaluate the effectiveness of the developed method in 
predicting the performance of statically operating propellers, calculations 
were performed for several propellers for which static data were available. 

The first two configurations, Propellers 1 and II, were tested at Texas A & M 
University and the measured data were reported in Reference 11. The third 
configuration. Propeller III, was tested at the Wright-Fatterson Air Force 
Base and the experimental results are reported in References 12 and 13. 

Reasonably good correlation between the thrust and power predictions 
of the new method used and the measured data resulted in the case of 
Propellers I and II. Correlation for Propeller III was not so good as the 
other cases. However, according to Borst and Ladden 14 , the Wright 
Patterson whirl rig that was used has a large cross-sectional area for the 
test rig relative to. the area of the propeller disk. This kind of situation 
would certainly influence the wake geometry and could very well cause 
inaccuracies in the prediction method, if the wake geometry built into the 
method were used. Although an attempt wa/5 made to alter the description of 
the wake geometry by allowing for blockage, no great difference in the 
calculated results was realized and it was concluded that it would be 
necessary to generate empirical wake data for this particular case. 

As a final part of -the work with this prescribed-wake method, a sensi- 
tivity study was made relative to the various parameters used in defining 
the geometry of the ring stacks. It was found that the predicted perfor- 
mance is most sensitive to the axial spacing ox the vortex rings. It was 
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found that the rate of contraction of the vortex rings had a relative effect 
in the thrust and power coefficient values but not in the Figure of Merit 
values. Although ways of improving the method were suggested in Reference 8, 
generally the computer program was found to be satisfactory for its intended 
purpose. 

Unsteady Vortex Lattice Technique Applied to the Wake Formation 

As previously noted, the principal record of the work with the unsteady 
vortex lattice technique was in Reference 7. In his work, Hall treated the 
blades aS lifting surfaces, in fact, the treatment was general enough to 
handle either a‘ propeller blade or a finite aspect ratio wing with only 
slight changes in the computer program. 

The numerical model of the lifting surface and its wake consisted of 
replacing the continuous distribution of vorticity by a mesh of vortex 
segments of finite length and strength. The geometry of the wake vortices 
was fixed by the motion of an ever increasing number of points moving under 
the influence of the bound vorticity and its own self-induced effect since 
it was assumed that these wake points are connected by straight-line vortex 
segments Identified as shed and trailing vorticity. The description of the 
blade bound vortices was fixed by the blade geometry. 

The vortices on the surface were arranged in a conventional manner. 

The surface was broken into a number of spanwise segments and each spanwise 
segment was subdivided into a number of chordwise segments. ‘ Each resulting 
panel contained a control point and was spanned by a straight-line vortex 
segment. The spanning vortex was at the 1/4-segment chord of each panel and 
the control point was at mid-segment span and 3/4-segment chord of each 
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panel. The final spanwise filament was 1/4-segment chord downstream of the 
surface trailing edge, implying that the Kutta condition was satisfied 
approximately, the accuracy of approximation increasing as the number of 
chordwlse segments increase. 

The surface loading was reflected in the strength of the bound vortices 
and the unique solution to the load distribution was determined by applying 
the boundary conditions of tangent flow at the surface and the Kutta condi- 
tion. The solution was obtained numerically by expressing these conditions 
as a system of simultaneous algebraic equations and solving by matrix multi- 
plication methods. The velocity associated with the vortices was described 
by the Biot-Savart law with the load distribution on the surface being 
determined from the unsteady Bernoulli equation. 

At the start of this study considerable time was spent with the finite 
„wing since this configuration contains essentially the same numerical 
.problems as the propeller but is .less complex. The first configuration 
_ithat was tried was the infinite wing. This was done by making the aspect 
ratio sufficiently large (AR ■ 1000) that it adequately represented the 
infinite aspect ratio case. The comparison with the Wagner solution was 
quite good except in the initial instants where large deviations occur. 

The explanation for this lies in the fact that the Wagner solution contains 
only the effect of the wake whereas the numerical solution contains an 
"Infinite" added mass ll't solution with the impulsive start. 

With regard to the finite wing, such things as the effect of the number 
of spanwise panels on the lift coefficient were investigated. Linearized 
wakes and wakes which distorted under the influence of velocities induced 
by the shed, trailing, and bound vortex systems were also studied. It was 
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found generally that, for the distorted wake case, there was little tendency 
for the vortices to roll up into tip vortices unless localized Induction 
was taken into account. The extreme slowness of the roll up rates was 
remedied by the use of localized induction concept 1 " 1 which, in effect, says 
that a curved vortex filament induces at a point on itself a velocity propor- 
tional to the local curvature and is directed along the local binormal. 

This means that the curvature of the trailing vorticity can induce a span- 
wise flow which will tend to destroy the initial two-dimensional character 
of the motion. Under this influence the vortex segment end points describing 
the wake will travel spiral paths which promote interference between 
filaments and increase the roll up lath. It was concluded that this 
localized induction effect is an essential ingredient as iar as a realistic 
roll up process in concerned. 

The Biot-Savart equation contains a singularity, if the point at which 
the induced velocity is determined lies on the filament. In order to circum- 
vent this, the common practice is. to assume a core exists in which the fluid 
moves as a solid body and not *as potential flow. Studies were made by Hall 

t t- 

with regard to the proper core size to use. Along with this kind of study 
Hall also assumed that the circulation about a vortex segment varied with 
the length of the segment as it distorted. It is with the idea of core 
size and the circulation as the segments distort that Daso 9 was primarily 
interested. 

With the analysis established and verified for the finite wing, Hall 
undertook the case of the statically thrusting propeller. He tried 
predictions of a four-bladed propeller whose performance was presented in 
Reference 16. It was found that the theoretical results for the propeller 
configuration did not correlate well with the experimental results. In an 
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effort to obtain further comparisons, a classical Prandtl analysis was 
performed and calculations based on momentum theory were made. In general, 
reasonable comparisons in thrust predictions were obtained between Hall's 
analysis and the Prandtl analysis while the actual measurements of 
Reference 16 were considerably lower. Then, using an average thrust 
coefficient, C^, of the value predicted by either the present analysis or 
the Prandtl analysis, a momentum power coefficient, Cp^, was calculated. 

It was found that Hall'v. nalysls compared favorably with this Cp^ value 
as well as with the total C p of the Prandtl analysis. However, all of 
these calculated values are much higher than the measured C p of Reference 15, 
indicating that possibly the measured C p was too low. The error observed 
in these results was much greater than anticipated, particularly since the 
finite wing results were so encouraging. 

Further error in the analysis could have been due to poor synthesis 
of the airfoil section data. Although care was taken and the guidelines 
of Reference 16 were followed, the airfoil section was nonstandard and 
difficult to describe. Poor estimates of the drag characteristics could 
explain, in part, discrepancies in the power calculations among analyses 
with reasonable thrust comparisons. 

Error might also have been due to the relatively short wakes generated. 
Even though extremely long computational run times (20,000 sec.) were 
performed, only about two revolutions of wake could be generated at best, 
and it is quite conceivable that this is not enough to predict the steady 
state performance. It was noted that the average C^, and Cp^ responded to 
the impulsive start much like a low aspect ratio wing. That is, following 
the impulsive start the performance dropped very quickly to what appears to 
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be the steady state value. It is possible that steady state had not been 
attained and more revolutions were necessary. This would lead to an 
increased inflow which, by decreasing blade angle of attack, could lead to 
decreased thrust prediction. Regions of inboard stall would be determined 
by this inflow, and performance would be measurably affected by the extent 
of these regions. 

Finally, there is an error due to the vortex wakes deposited by the 
propeller blades. The time steps considered were generally much too large 
to predict accurate wakes. As a result, the wakes of the four-bla.ded 
configuration became unstable; this instability was enhanced by interaction 
core radii that were too small. The resulting wake geometries then contained 
extremely long straight line vortex segments which, once formed by a strong 
interaction induced velocity acting over a relatively large time step, 
could produce completely erroneous velocities at the blades. To make matters 
worse these segments could never return to a reasonable geometry as time 
progressed since they might never pass through" enough interactions to 
counteract the effect of one strong one. It should be noted that wake 
instabilities noted in the analysis are believed to be only numerical with 
no physical counterpart. 

Even though the comparison of theoretical and experimental results 
leaves much to be desired, some parametric results were successfully obtained. 
Small time steps (1.5° to 3° in azimuth) are required for accurate wake 
prediction. This is necessary to determine an accurate vortex filament radius 
of curvature for calculating the locally induced effects. This is also a 
requirement in order to obtain reasonable vortex induced curved paths for 
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the wake points from the one-step Euler integration scheme which can only 
provide straight line translation of a point. 

One of the important regions of concern is the location where the 
wake from a preceding blade comes close to the following blade. .The accuracy 
of these blade-wake interactions not only depends on small time stepB, bu.. 
also on interaction core radii large enough to limit the movement of a wake 
point to a reasonable value. 

The conclusion of this work must admit that the accuracy of the present 
analysis when applied to the statically thrusting propeller has not been 
satisfactorily demonstrated since correlation with the selected experimental 
results was poor. Even though the basic formulation is believed sound from 
comparison with other analyse^ and finite wing results, final correlation 
will have to await better experimental results, more accurate airfoil section 
characteristics, relief from the numerical inaccuracies associated with the 
aerodynamic interference region and larger computational runs to numerically 
establish the wake. This procedure, like other vortex lattice techniques, 
uses an Inordinate amount of computer,’ time due to the repeated calculations 
of the Biot-Savart Law in the wake. Unfortunately, no wake simplification 
or approximations are apparent because of the importance of the nonlinear 
flow of the induced velocity field. This is further aggravated by the 
small time step requirement to compute interference aerodymamlcs of the problem 
accurately. This seriously restricts the usefulness of the analysis, at 
present even as a research tool. However, vortex lattice techniques are 
those which most readily apply to nonlinear aerodynamic problems so that 
further attempts at reducing the computation time of this analysis, as 
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well as accepting long time computer runs, are perhaps justified, at least 
in reseurch problems*, 

In spite of thy inconelusiveness of the primary results of this 
analysis, some positive results were obtained. Perhaps the most significant 
of these is the modeling of the wake roll~up with the localised induction 
concept while considering the three-dimensional flow about a lifting 
surface starting from rest. 

Vortex Core Size Study 

As mentioned earlier, the vortex core size and how the circulation 
varies with elongation of the vortex segments was of concern in this general 
study. Daso^ was concerned with providing some rational approach to 
determining core sizes which wasn't just an arbitrary choice. He was also 
concerned with keeping track of core sizes as the segments distorted and 
this aspect of the problem was intimately related to what happens with 
the circulation during distortion. 

In general he showed that the circulation must remain constant regard- 
less of segment length. On the other hand, because of the conservation of 
mass in the core, the vorticity will increase with increasing length of the 
vortex segment and the core radius will correspondingly decrease. Therefore, 
once having established a core size, it is a matter of routine to keep 
track of the core size as the vortex segments distort. 

With regard to establishing the initial core radius, Daso first looked 
into an approach which took cognizance of the fact that at the trailing 
vortex sheet the velocity induced just above and just below the sheet are 
proportional to the circulation per unit length. By dividing the sheet 
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Into sections, the circulation per unit length is known from the circula- 
tion distribution and in turn the induced velocity is known. The vortex 
core size can then be found via the Biot-Savart law, but since the induced 
velocity depends on the original arbitrary chord of sheet segment size, 
the vortex core that results is also arbitrary. 

Daso then studied the pressure distribution for a Rankine vortex with 
the thought that, by using reasonable values of minimum pressure at the 
center of vortex, the core radius could be explicitly calculated. Although 
the minimum pressure at the center of the core is generally unknown and 
cannot be theoretically determined, plots of core radius versus minimum 
pressure show a range of pressure where little change in core radius occurs. 
Beyond this range the core radius increases quite rapidly and a zero 
minimum pred<‘uf>e in the other direction is an improbability. On this basis 
an average value of minimum pressure of about 400 lbs per square foot was 
chosen. The error in radius involved in this choice varies from about 
8 percent at a minimum pressure of 100 lbs per square foot to 15 percent 
at 800 lbs per square foot. Although it was recognized that there was 
a degree of arbitrariness in the choice of core size, it was reasoned that 
this degree was relatively small. With the minimum pressure chosen, the 
core size becomes a function of circulation. Daso made studies of how the 
initial estimate of core radius varied with forward velocity, spanwise 
position, and with time step. He did this for both the trailing and shed 
vortex segments". He did, however, restrict himself to the use of the 
program with wings only because, as in Hall's work, it is much easier to 
work with wings rather than the propeller when proofing such features as 
just discussed. 
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The question of how the changes made in the program would affect the 
final results when running the propeller case is still unanswered. Time 
and money did not permit such an exercise. A listing of the general 
computer program is presented in the Appendix. The code contains comment 
statements for parts of the program that were used in the vortex core study 
just described. These can be included in the program by just eliminating 
the comment designation. The same can be said for a number of other . 
Statements which give the option of running program with the IBM 370 
facilities at Penn State or at the NASA Langley facilities. 

Conclusions 

In spite of the generally favorable trends established from applying 
vortex lattice techniques to the statically thrusting propeller, the 
primary objective of obtaining the high degree of accuracy necessary to 
correlate theory and experiment has not been accomplished. However, 
the major problem areas in the aerodynamic modeling have been identified 
and the foregoing analysis represents a tool to investigate these areas. 

If more frultfull results are to be obtained, efforts to reduce 
computer time must be of prime concern. Attempts to more accurately 
predict the potential Inflow lead to small time Increments corresponding 
to an azimuth step size, A0 < 1.5°, fully one-half the smallest value 
considered and at least one tenth a value at present practical. This 
limit has beer, established by estimates necessary to promote good wake 
roll-up characteristics. Attempts in the present analysis to reduce 
computer central processor time (and vore storage) with special data 
handling techniques have been generally unfruitful. 
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Reductions in computation time would also permit more accurate 
representations of the wake. The numerical integration scheme considered 
in the present analysis is a simple one-step Euler scheme, shown to be 
less accurate than either a Runge-Kutta method or a one-step predictor- 
corrector technique. The inherent inaccuracy of the method lies in the 
fact that points can only translate under the influence of a vortex 
Induced velocity whereas the true path is circular. Unfortunately this 
method is the most economical from the point of view of computation time 
and core storage, although to get a sufficiently close approximation to 
the circular path requires very small time increments. 

The final work on core size and the variation of vortlcity in the 
core appears to have resulted in a satisfactory means of handling these 
quantities. 
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APPENDIX 

Computer Program for the Unsteady Vortex Lattice Technique 

The following is a listing of the computer code used in the unsteady 
vortex lattice approach. As can be seen in the following description 
of input data, certain choices permit the program to be run for a wing 
or a propeller. For instance, the choice of zero rpm and one blade 
permits the analysis of a wing. 

First Data Card 

NOPAN - Number of spanwise panels of the lifting surface. 

NUM - Number of spanning vortices including the shed vortex on 
each spanwise panel. 

MXTIME - Maximum number of time steps. 

IBL - Number of blades. 


Second Data Card 


V - Forward or free-stream velocity. 

RPM - Revolutions per minute. 

R - Radius of blade or wing span. 

BL - Number of blades. 

DELTH - Angular increment of blade travel in one time step. 
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Third Data Cards 

YC(L) - Control point coordinate along spanwise or Y-axis of Lth 
control point. : 4 

BETAC(L) - Pitch angle at a control point or angle of attack of wing 
at the Lth control point. \ 
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THC(L) - Section thickness at a control point. 

CC(L) - Chord length of: a lifting surface element corresponding to 
Lth control point. 

XRC(L) - X location with respect to the leading edge of the stackup 
point or origin of the blade based coordinate system. 

Fourth Data Cards 

The input variables, Y(L), BETA(L) , TH(L) , C(L), and XR(L) have 
similar definitions to those of the third data cards except that they 
refer to the edge of the lifting surface panel containing the control 
point, L. 

Other Incut Variables 

i — -■ 1 1 - - » i«cas- 

RHO - Nondlmensional fluid density. 

ROR - Density of air, 

A sample of the input data used in the case of a wing study appears 
on the last page, following the program listing. The occasional 
printing of CANADAIR in the program is in -'oference to propeller section 
data for a CANADAIR propeller. 
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C PROGRAM PENNSTI INPUT, OUTPUT, TAPE5= I NPUT , TAPE1 ) 

DOUBLE PRECISION DARG 
C REAL IT I MUZ 

INTEGER STATUS 

DIMENSION 0 { 100 , 1 ) F IPIVOT( 100 ) , INDEX ( 100 , 2 ) 

DIMENSION A( 100, 100 J r F( 100) 

DIMENSION XC< 10,30 ) , YCI30) -,B( 120) , VN ( 120) , TGAM ( 120 ) 

DIMENSION ZC{ 10,30 ) , WXC ( 1C , 30 ) , WYC l 10 , 30 ) , WZC ( 10 , 30 ) , FXQS ( 10 , 30 ) , 

1 FYQStlO, 30), FZGS( 10,30) , FXUS { 10 , 30 ) , FYUS ( 10 ,30 ) ,FZUS ( 10 ,30 ) , 

2 THRST ( 30 ) , DRAG I 30 ) » TORQ ( 30 ) 

DIMENSION XK( 30) ,XRC( 30 ) ,DZP ( 20,30) , TH( 30) ,THC(30) 


DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
C CM MON 
100 ) , GAMMA { 100 ) , GAM 
10) , BETA ( 30 ) , AS I 30) 


VI (30,50 ) , VJ I 30,50) , VKI3 0, 50) 

FXBV ( 10 , 30 ) » FYBV { 10,30),FZBV(10,30) 

PQS( 10,30) , PUS f 10,30) ,P( 10,30) ,P0WR(30) 

PM IN ( 30) » SUMAR ( 30 ) 

D { 3 ) » STA TUS ( 4 ) 

S I NLAM (10,30) , COSL AMI 10,30 ) 

X( 10, 30) , Y( 30), 2(10,30), XW(30, 100) , YW ( 30 , 100 ) , Z W ( 30 , 1 
S { 30 , 100 ) ,GAMT(30,100) , AN ( 30 ) ,RA ( 30 ) ,0MA(3O) ,C(3 
,U( 30,31) » YYC ( 30 ) ,CC(30) ,BETAC(30) »RAS(30,100) , 


1000 

1001 


1RAT( 30, 100) 

COMMON VX P » VYP , VZP , CO T , S I T , I T IME , I WAKE » BL , I BL , NOPAN , NUM , XD , YD , ZD , 
1H, E,AO,R,KMAX,LlNWA,V,SUMARR 
DATA STATUS/4*0/ 

FORMAT ( 1H ,9( E13.5 ) ) 


FORMAT ( 1H , 5X, 'THRUST DISTRIBUTION' 
1' EFFECTIVE ALFA DISTRIBUTION') 


10X, 'POWER DISTRIBUTION', 7X, 


1002 FORMAT ( 1H 

1 ' PROPELLER 

1003 FORMAT ( ' ' , 
FORMAT ( 1H0 ) 
FORMAT ( ' 


,10X, 'THRUST COEFFICIENT' , 10X, 'POWER COEFFICIENT' , 10X' 
CONVENT ION' ) 

i 


• , E13. 5 ) 


1004 

1005 

1006 

1007 

1008 

1009 

1010 
1011 
1012 

1013 

1014 
1016 
1017 


1 5 , 6 { El 3 • 5 ) , 15) 


110 ) 


FORMAT ( ' 0 ' , 9 ( E 13 .5 , IX ) ) 

FORMAT ( »0» ,2( 15, 2X ) ) 

FORMAT ( ' ' , 100X , *T IME USED=' 

FORMAT ( ' ' ,5( 10X,E 15.8) ) 

FORMAT ( 'O', 10X » 2 ( E 13 • 5 ,5X ) ) 

FORMAT ( 'O', 10X , E 13 . 5 ) 

FORMAT (715) 

FORMAT ( 7F10.6 ) 

FORMAT ( 1 MAXIMUM NUMBER OF TIME STEPS IS 

FORMAT ( 'O', L 2 { ’CAN ADAIR M) 

FORMATJ *0* , 'NUMBER OF SPANWISE P ANELS= ' ,13//' 


VORT I CES= ' 


IE 

2 BLADE NUMBER** 1 , ’12). 

1018 FORMAT ( 'O' , 'FLIGHT SPEED=' 


13//' 

f , 


MAXIMUM NUMBER OF TIME 


,15) 

NUMBER OF CHGRDWISPR0P0465 
STEPS=* , 


14// 


E12.5, 'FEET PER 


1E12.5//' RADIUS! SPAN ) = ' ,E12. 5, ' NUMBER 
2 DELTA THETA=* ,E 16.8, 'DEGREES' ) 

TIME I NCR EM ENT = ' , E 1 6 .8// ' 


SECOND R PM= * 

□F BLADES= ' , E12. 5 // ' 


1019 FORMAT ( ' 0 • , ' 
1PANEL=* , £12.5/7 ' 

1020 FORMAT ( 'O' , • 


1 //' 

2 //' 


YC/R 


REL. LENGTH 
COEFFICIENT-' ,E16.8) 

BLADE SECTION CHARACTERISTICS 
CONTROL POINT GEOMETRY 
BETAC ( DEC . ) TC/C 


OF CHORDWISE 


PR0P0470 
PR0P047 1 
, PR0P0475 
PRDP0480 
PR0P0485 
PR0P0490 
PRPP0495 
' PR0P0500 
' PR0P050 5 
PR0P0506 


ORIGINAL E IS 
OF POOR OVALITY 


oooooooonoo 
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3CC/R 

1021 FORMAT ( 

1022 FORMAT ( 
l//* 

2C/R 

1023 FORMAT ( '0 
1050 FORMAT ( ' 
7777 FORMAT ! 1H 


•O' , 

Y/R 


' ,6(E13 

, 


T/C 


XRC/R* ) 

5 , 9X ) ) 

VORTEX GEOMETRY 
BETA ( DEG. ) 

XR/R* ) 

,' BOUND VORTICITY DISTRIBUTION'/) 
i •*** PH I = * * E 16 . 8 ) 

,'SHEO VORTEX CORE RAD I US , RAS ' , 5X , ' TRA I L I NG VORTEX COR 


PROP0507 

A 


' PR0P0515 
PROP05 16> 
PR0P0517] 
PRDP0520 ; 




5555 

6666 

8888 


IE RADIUS, RAT' ) 


FORMAT! 'O' , 'MINIMUM PRESSURE, PMIN LBS PER FT**2') 


FORMAT! ' 0 », 20X, 2! E 13.5,20X1) 

FORMAT! 1H ,3! E13.5 , 20X ) ) 

LINEARIZED BOUNDARY CONDITION REMOVED 
LINEARIZED OR DEFORMED WAKE 

LINWA-1 IMPLIES LINEARIZED WAKE 
LINWA=ANYTHING ELSE IMPLIES DEFORMED 


WAKE 


NOTE** IN FREE WAKE ANALYSIS, CONSERVATION OF CIRCULATION HAS 
GAMT*AL=CONST AT SHEDDING !FUNC. ONLY OF TIME OF SHEDDING AND 
SPANWISW POSITION). 2! B-S LAW GAM=(GAMT*AL)/AL(T)=CONST/AL 
WHICH LEADS TO AL**2=ALS APPEARING IN DENOMINATOR INSTEAD OF AL 


CALL TIMUSE! ITIMUZ ) 

PRINT1008 , ITIMUZ 

C CALL BLKREWD(5LTAPE1) 

PRINT 1016 
LINWA=0 
C L INWA= 1 

READ I 5. 1012) NOPAN ,NUM,MXTIME, IBL 
PRINT 1017, NOPAN, NUM,MXTIME, IBL 
READ (5,1013) V,RPM ,R,BL,DELTH 
PRINT 1013»V,RPM,R,BL,DELTH 
RHO = 1 . 

BL= I BL 

PI=3. 1415927 

C**** MDIM MUST BE GREATER THAN OR EQUAL TO FIRST SUBSCRIPT OF ARRAY A 
C*** (MAIN PROGRAM), SO THAT ARRAY A (SUBROUTINE MXINV) IS PROPER 
MD I M= 100 
CED ITIMG=0 
HH=0 .001 
E=. 000001 
H=HH 

C**** V=0 IMPLIES HOVER 
C *** THET A=0 IMPLIES PSI =90 
THETA=0 . 

RPS=RPM*P 1/30. 

VT I P=RPS*R 
DELTH=DELTH*P 1/180 . 

COEF=RHO*PI*R*R*VT IP*VTIP 
IF(COEF.NE.O.O) GO TO 107 
CCEF=0.5*RH0*V*V*R*R/3.0 
C DELT = D5LTH/R,PS 


li 

4 

107 DELT=0.1 

C**#************nOPAN-NUMBER df spanwise panels 

C******#******** NU M:=NUMB ER OF CHORDWISE VOR T I C E S , I NCLUD I NG SHED AT 
C*#******#*#*** H XTIME=MAXIKUM NUMBER DF TIME STEPS IN WAKE 
Cfc***^-***********#V = FLIGHT SPEED 
C?*»r*#**$ c******** RPM=R PM 
C **»!>*** *#*#*#**«*# R=RAO I US T SPAN 
Css*************** BLr IBL=NUMBER DF BLADES 

c******* delth=angluak increment df blade travel in dne time step 

C **********************£ hq=flui d density 

c********* coep=nondimensionalization factor for forces 

q******************** OELT s TIME STEP INCREMENT ******************* 
NUMM I=NUM-1 
NPANP 1=N0PAN+1 
NPANM1=N0PAN-I 
NPANPA = NOPAN + 2 
MATRX1=NUM*N0PAN 
MATRX3=NUMMI*N0PAN 
MATRX2=MATRX3+ 1 
A0=0. 1076 
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i ■ 

i ■! 
! ; 


CANADA IS 


A 


DELX=1./NUMM1 

PRINT I0l9t DELT » DELX» COEF 

PRINT 1020 

DO 113 L=1?N0PAN 

READ ( 3 » 1013 ) YC( L ) » BETAC ( L ) , THC ( L ) , CC ( L ) ,XRC(L) 

PRINT 1021 f YC { L ) ,B ETAC ( L ) v THC ( L) tCC { L ) ,XRC (L) 

BETAC (L)=BETAC(L)*P 1/180. 

YC(L)=YC(L)*R 
CC(L)=CC(L)*R 
113 CONTINUE 
PRINT 1022 
DO 1 1 A L~lt NPANP1 

READ (S ? 1013 )Y(L) ,B ETA ( L ) , TH ( L ) ,C(L) ,XR(L) 

PRINT 1021»Y(L) »BETA(L),TH(L) »C(LJ >XR(L) 

BETA ( L ) = BETA(L)*PI / 1 8 0 . 

Y ( L ) =Y ( L ) *R 
1 1 A C(L) = C(L)*R 

DO 115 L=l, NOPAN 
DO 115 1= I? NUMM 1 

115 DZPU,L) = 0. 

C ************** :*****s*nxCfWYC T WZC INITIAL I ZED*************************** 

DO 116 L=l, NOPAN 
DO 116 1 = 1 1 NUMM 1 

WXC( I rL ) =0 . 

WY C ( I r L ) =0 . 

C 100 

116 WZ C ( I , L ) =0 . 

C*** XC IMPLIES CONTROL POINTr VORTEX COORDINATE WRT BLADE SYSTEM 
DO 3 L=l, NOPAN 
DO 3 I = 1 » NUMM1 
R I = I 

XC( I, L )= ( (RI-.2S )*OELX-XRC( L ) )*CC ( L ) *COS ( BETAC ( L ) ) 

3 ZC( If L)=-( (RI-.25) *OELX-XRC(L ) )*CC(L)*SIN(BETAC(L) ) 
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DO 2 L«1»NPANP1 
DO 2 I=1,NUM 
R 1= I 

X ( I , L > = t (RI-.75I*DELX-XR(L) )*C ( L ) *CQS ( BETA ( L ) ) 

2 Z( I,L)=-t (RI-.75)*DELX-XR{L) ) *C ( L ) *S I N ( BETA ( L ) ) 

CODETERMINATION OF LOCAL DIHEDRAL**#* 

DO 600 1=1, NOPAN 
DO 600 I = 1 , NUMM 1 

ALAM=-ATAN( -XC ( I,L )*( TAN ( BETA ( L+ 1 ) ) -TAN ( BETA ( L ) ) 1 / ( Y i L+ l ) -Y ( L ) ) J 
C * ALAM=G 

SINLAMI I , L )=S INI AL AM ) 

600 COSLAM ( I,U-COS( ALAM) 

21 11=0 

C** COEFFICIENTS DETERMINED ON BASIS OF STRIP THEORY-NO SPANW1SE EFFECTS 
DO A L“l, NOPAN 
CCOSBL = CQS(BETA(L) )*CC(L) 

CSINBL=SIN(BETA(L) )*CC(L) 

DO A I = l, NUMM 1 
XCI = XC( 1 , L ) 

ZC I=ZC ( I , L ) 

YC 1= YC ( L ) 

XD=XC I 
YD=YC I 
ZD=ZC I 

C** ZN,XN,YN UNIT NORMAL COMPONENTS WRT BLADE FIXED AXIS 
ZN=COS( BETACILJ-DZPI I ,L ) )*COSLAM( I , L ) 

XN=SIN(BETAC(L)-DZP( I,L) ) *COSLAM ( I , L ) 

YN=S INLAMI I , L ) 

I 1=11 + 1 
J J = 0 

DO 4 K=l, NOPAN 
D04 J= 1 , NUM 
JJ=JJ+1 
A 1=0 . 

A2=0. 

A3 = 0. 

DO 1IB=1,IBL 

CSIN=C0S(2.*PI*( IB-D/BL) 

SSIN=SIN(2.*PI*( IB-D/BL ) 

DO 1 KKK = 1 » 3 
IF(J-NUM) 111,112,111 

111 I F ( KKK-2 ) 109,112, 110 
109 XB=X ( J,K)*CSIN-Y(K )*SSIN 

XA=X(NUM,K)*CSIN-Y (K)*SSIN 
YB=Y(K)*CSIN+X( J,K )*SSIN 
YA=Y(K)*CSIN+X(NUM ,K)*SSIN 
ZB = Z C J,K) 

ZA=Z ( NUM , K ) 

GO TO 106 

112 XB=X(J,K+1) *CS IN~Y (K+1)*SS1N 
XA=X( J,K)*CSIN-Y(K )*SSIN 
YB=Y(K+1)*CSIN+ X( J ,K+1 )*SSIN 
YA=Y ( K ) *CS IN+X I J »K 1#SS IN 


oonooooooooo 


ZB=2( J,K+1) 

ZA=Z( J f K) 

GO TG 108 

110 XB=X<NUM,K+l)*CSlN~Y(K+l)*SSIN 
XA=X(J»K+i)*CS IN-Y (K+D+SSIN 
YB=Y(K+1 )*CSIN+X(NUM,K+i )*SSIN 
YA=Y (K + l ) *CS 1 N + X { J »K+1)*SSIN 
ZB= Z ( NUM r K+ l ) 

ZA=Z ( J, K+L) 

108 CONTINUE 
XBA=XB-XA 
YBA=YB-YA 
ZBA=Z6-ZA 
XDA=XD-XA 
YD A=YD-YA 
ZDA=ZD-Z A 
XDB=XD-XB 
YDB=YD-YB 
ZDB=ZD-ZB 

ALS=XBA*XBA+YBA*YB A + ZBA^ZBA 
ACS=XOA*XDA+YDA*YD A+ZDA*ZDA 
BCS = XDB*XDB+YDB*YD B+ZDB*ZDB 
AL= SORT ( ALS ) 

AC= SORT ( ACS ) 

BC= SGRT(BCS) 

CO$A= ( AC5+ALS-BCS) /(AC*AL*2. ) 

COS B= ( BCS+ALS-ACS) /( BC*AL*2. ) 
TEMPA=l.-COSA*COSA 

117 VFN=(C0SA+C0SB)/(AL*ACS*TEMPA*Pl*4. ) 
AILXAC=YBA*ZDA-ZBA*YDA 
A JLX AC=ZBA*XDA-XBA *ZDA 
AKLXAC=XBA*YDA-YBA*XDA 
A1=A1+VFN*AILXAC 
A2=A2+VFN*AJLXAC 
A3=A3+VFN*AKLXAC 
IF(J-NUM) 1 » 4 » 1 

1 CONTINUE 

4 A ( I I » J J )=A1*XN+A2* YN+A34ZN 
GO TO 1193 
1193 PO = 2116.8 

ROH = 0.002378 

PM I N = 600.0 

00 1191 Ji=l»10 

PM IN I J1 )«J 1*100.0- 100.0 

SUMAR { J 1 ) =SQRT( (PD -PM INI J 1 ) )/ROH) 

SUMARR=SUMAR ( J 1 ) 

SUMARR = SORT ( ( PO - PMIN ) /ROH ) 
PRINT 6666 » PM I N ( J1 } 

DO 88 L= 1 » NPANP 1 

NPANP3=L 

SUMD=0.0 

DO 87 N= l » NP ANP3 

BA=2.0*N-1.0 


o r> o non rtoooo 
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SN=SQRT ( BA ) 

87 SUMD= SUMD + AN ( N ) *SN 
RAIL) a ( V*R *SUMD)/SUMAR(J) 

PRINT 8008, AN ( L ) , R A ( L ) 

88 CONTINUE 
KKK=1 

DOll I 1=MATRX2, MAT RX L 
LLL=KKK+NUH 
D012 J J= 1 * HATRX l 
All I , J J ) =0. 

IFUJ.GE.KKK.AND.JJ.lt. ILL) A(II,JJ)=l. 

12 CONTINUE 
il KKK=KKK+NUM 

DO 100 I=1,MATRX1 

100 PRINT 1000, ( A ( I,J) , J=1,MATRXL) 

CALL MAT INV ( A, MATRX1,Q,0, DETERM, IPIVOT, INDEX, 100, I SCALE) 

CALL MX1NV(A,MD1M* MATRX1) 

PRINT 1011, DETERM 
PRINT 10 12, ISCALE 
E PMIN=0.0 

PMIN=400.0 
PRINT 5555 
PRINT 6666, PMIN 
C3003 CONTINUE 
I T I ME=0 

SUMARR=SQRTl ( PO-PM IN) /ROM) 

71 CONTINUE 

CALL T IKUSEI IT I MUZ ) 

PRINT1008, ITIMUZ 

RTIME=ITIME 

TIME=RTIME*DELT 

THET A=D ELTH*RT IME 
SINTH=SIN(THETA ) 

COSTH=COSlTHETA) 

C********** IMPACT VELOCITY — ITIME-O ********* 

C *** COMPUTATION OF IMPACT VELOC I TY , VN I 1 1 ) *** 

11 = 0 

DO 25 L=l, NOPAN 
DO 25 1=1, NUMM 1 
I I^TI+I 

25 VN( II )=(RPS*YC(L)+V*CGSTH)*SIN<BETAC(L)-DZP(I,L) ) *COSLAM ( I , L ) + ( -RP 
1S*XC( I,L)-V*SINTH) *SINLAM( I,L) 

IF! IT IME ) 14, 19,20 

19 D06 11=1, MATRX3 
6 F ( I I*)=-VNC 1 1 ) 

DO 27 II=MATRX2,MATRX1 
27 F ( 1 1 ) =0. 

GO TO 29 

20 CONTINUE 

C**** WAKE BOUNDARY VELOCITIES WRT BLADE-FIXED SYSTEM ************** 

11 = 0 . 

PRINT 7777 
D08 L=1,N0PA'! 


CED 


<~»o n n nonnnorto 
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DOB I=1,NUMM1 ! 

XD=XC(I,L) j 

YD= YC ( L ) I 

ZD=ZC l I » L ) ! 

11 = 11+1 | 

C****** back transform wake coordinates to blade -fixed system ************** 

COT-COSTH 

SIT=-SINTH 

I wake~ 1 

H = HH 

CALL INVEL 

I F ( KMAX • EQ.NUM ) GO TO 1149 
I F ( Kf'.AX . EQ.NPANP l ) GO TO 1149 
IF(L.GT.l) GO TO 1 149 
IFU.GT.l) GO TO 1149 
DO 1132 Kl=WNPANPl 

1132 PRINT 86S8> AN(Kl), RA(K1 ) I 

1132 PRINT 8bOB , KAS ( K1 # IT IME+i ) , RAT { K 1 » IT I ME+l ) 

PRINT1006»VXP,VYP» VZP j 

1149 WXC( I * L ) = VXP 

wyc( ifLi-vYP 

WZC( I »L)=VZ.P ii 

**** BOUNDARY CONDITION FROM STRIP THEORY j 

8 FI I I )=-VN( II )-VZP*COS(BETAC(L )-DZP( I ,L) ) *CQSLAM ( I , L ) -VXP*S I N ( BET AC 
1 { L ) -D2P C I,L) )*CQSLAM( I , L ) -VYP*S I NLAM ( I,L> 

************************ kUTTA CONDITION ******************** 

I I = MATRX3 - 

D016 1=1, NOPAN % 

11=11+1 I 

M=(L-1 )*NUM+1 ; 


N=M+NUM-2 
SUM=0« 

DO 17 J J=M, N 

17 SUM=SUM+GAMMA(JJ) 

16 F ( I I ) =SUM 
29 CONTINUE 

PRINT 1004 

PRINT1000, IFC I ) , 1= 1 , MATRX 1 ) 

CALL MXMLTI A, F, GAM MA , MATRX1 , MATRX 1 , 1 , 100 , 190 , LOO ) 

C **************** BOUND VORTICITY OUTPUT ***************** 
PRINT 1023 
DO 18 L=l, NOPAN 
M= ( L — 1 ) *NUM+l 
N=M+NUM-1 

18 PRINT 1000, (GAMMA! 1 1 ) , I I = M,N) 

PRINT 1004 

C******FORCE DUE TO DELTA P-CHDRDWISE LOADING 
C ****** FORC E ^ DUE TO DELTA P-CHORDWISE LOADING 
C**** FORCE ON PANEL DUE TO DELTA-P 
DO 33 L= 1 , NOPAN 
M= ! L- 1 1 *NUM+ 1 
N=M+NUM-2 
1=0 


4 


l 
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C*** M f N LOCATE CONTROL POINTS AND PANELS 
C*** LrI=CONTROL POINT INDICES 
DO 34 I I=M» N 
1 = 1 + 1 

C*#* DETERMINATION OP QUASI-STATIC FORCE ON PANEL ( I * L ) 

C*** 1 GAHI = CHORDWIS6 VORTICITY OT LEFT OF CONTROL POINT 
C *** GAM3=CHORDWISE VORTICITY TO RIGHT OF CONTROL POINT 
C*** GAM1 , GAM3 l + ) FEEDING INTO TRAILING EDGE 
GAM 1=0 . 

GAM3=0. 

C GO TO 601 

DO 137 J=M» 1 1 
IF(L-l) 138,139,138 

139 GAM1=GAM1+GAMMA( J) 

GAM3=GAM3+GAMMA( J+NUM) -GAMMA! J ) 

GO TO 137 

138 I F ( L-NOPAN 1 L4 1 , 140 , 14 1 

140 GAM1=GAM1+GAMMA ( J) -GAMMA ( J-NUM > 

GAM3 S GAM3-GAMMA ! J) 

GO TO 137 

141 GAM 1=GAM 1+GAMMA ( J ) -GAMMA ( J-NUM ) 

GAM3=GAM3+GAMMA( J+NUM )‘-GAMMA( J ) 

137 CONTINUE 
131 CONTINUE 

G AM2= GAMMA ( 1 1 ) 

C*** GAM2 ! + 1 LEFT TO RIGHT 

C ******************** INITIALIZATION OF FORCES ***************** 
FXQS ( I , L ) =0 • 

C 300 

FYGSU ,L )=0 . 

FZQS ( 1 1 L ) =0 . 

PQS ( I , L ) =0 - 
DO 142 KKK= 1,3 
IF(KKK-2) 144,145,146 

144 XD=(X( 2 ,L)+X(I+i,L ) ) / 2 . 

YD=Y ( L ) 

ZD=(Z( I ,L)+Z( I+1,L ) ) / 2 . 

J=I 
K = L 

JJ=I+1 

KK=L 

GAMA=GAM1 
GO TO 148 

145 XD= ( X ! I , L ) +X ( I , L +1 ) )/2. 

YD=(Y(L)+Y(L+1) )/2. 

ZD= { Z I I , L ) +Z ( I , L+l ) )/2. 

J = I 


146 


K = L + 1 
JJ = I 


KK = L 

GAMA=GAM2 
GO TO 148 

XD=(X( I,L+l)+X( 1+1 ,L+l) )/2. 
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YD — Y(L+l) 

ZD- ( Z i I ,L+1)+Z( f+l ,L+l) )/2. 
j = i 
K = L + i 
JJ=I+I 
K K = L + 1 
GAMA=GAM3 

148 I WAK E=0 
VX-KPS#YD+V*C05TH 
VY=-RPS*XD-V*SINTH 
VZ = 0. 

C***********VEL0C ITIfcS DUE TO BOUND VORTICIES CALCULATED IN BLADE -FIXED SYSTEM*, 
COT= I . j 

S I T = 0 . I 

H= .00000 1 1 

C***.* ************* ******************* **** ***** 

151 CONTINUE 

CALL INVEL \ 

I F ( KMAX . EQ.NUM ) GO TO 1150 
IF(KMAX.EO.NPANPl) GO TO 1150 
IF(L.GT.l) GO TO 1150 
I F ( I I.GT.l) GO TO 1150 
IF(KKK.GT.l) GO TO 1150 
PRINT 7777 
DO 1133 K2= 1 ? NPANP 1 

1133 PRINT 8888 T ANIK2)f RA(K2) \\ 

L 133 PRINT 8888, RAS ( K2 , I T I ME+ l ) , R AT { K 2 , I TI ME+ 1 ) < 

1150 VX=VX+VXP ;j 

V Y= V Y+VY P ! ! 

VZ=VZ+VZP 

VP0WX = 2 . *( RPS*YD+V*COSTH)-VX 

VP0WY=2 . *{ -RPS*XD-V*SINTH)-VY ') 

VPOWZ--VZ 

I F C ITIME.EQ.O) GO TO 150 

IF(IWAKE-1)149, 150 , 150 $ 

149 I WAKE=1 - 1 

C***** BACK-TRANSFORM WAKE COORDINATES TO BLADE-FIXED SYSTEM ******************* 

COT =COSTH . 

S I T = -S INTH 

H=HH ' 

C ************ ************** ************** i 

GO TO 151 •! 

150 CONTINUE 1 

FXQS ( I , L ) = RHC* ( VY* ( l ( J , K > - Z ( J J , KK ) ) - VZ* ( Y ( K ) -Y ( KK ) ) ) *GAMA . I 

1 + FXQS ( I » L ) -J 

FYOS ( I , L ) = RHO*( VZ* { X ( J , K ) -X { J J , KK ) ) -VX* ( Z ( J , K ) -Z ( J J ,KK )) ) *GAMA 
1+FYQS ( I »L ) . 

FZQS ( I »L J=RHO*( VX* (Y(K)-Y(KK) ) -VY* { X I J , K ) -X l J J , KK ) ) ) *GAMA 
1 + FZOS ( I ? L ) 

PQS ( I , L ) =F XQS ( I ,L )*VPOWX + FYOS { I » L } *VPOWY + FZQS { I ,L)*VPOWZ 
1 +PQS ( I , L ) 

C PRINT9999, I,L,KKK, XD, YD,ZD,VKfVY,VZ,GAM,A 

9999 f- 0 R M A T { 1 * , 3 ( IX , 13 ) , 7 C IX , E15 . 8 ) ) - 


14Z CONTINUE j.j 

C »*#«**♦ PRECEDING quasi-static forces have leading EDGE SUCTION *********** |j 
C*** DETERMINATION OF UNSTEADY FORCE ON PANEL! It L) If ' ij 

SUMP=0. | 

DO 35 J=M,II || 

.j IF! I T I HE ) 36,36,37 || 

36 SUMP=SUMP+GAMMA( J) J 

GO TO 35 M 

37 SUMP=$UMP+GAMMA( J) -TGAM! J) || 

35 CONTINUE f 

SUMP=SUMP/DELT*RHD | 

C*** UNIT NORMAL COMPONENTS AT XC(I,L) FROM ( AL ) X ( AC ) r ( AL*AC ) WITH I 

C*** AL=S PANWISE VORTEX SEGMENT AND AC=(L+1) CHORDWISE SEGMENT f 

XBA=X( I, L+D-Xl I , L ) t 

YBA=Y(L+1)-Y!L) | 

ZBA=Z! I ,L+1)-Z! I,L ) 

XDA=X( I » L+l )-X( l + l , L+l ) | 

ZDA=Z{ I , L+l )-Z ( 1+1 tL+1) 

C AL= SORT t XBA^XBA+Y BA*YBA+ZBA*ZBA) -- '[ 

DARG=XBA*XBA+YBA*YBA+ZBA*ZBA 
AL^DSORT ! DARG ) 

C AC= SQRT(XDA*XDA+ZDA*ZDA I 

DARG=XDA*XDA+ZDA*ZDA 

AC=DSORT { DARG ) j i 

AILXAC=YB4*ZDA j 

AJLXAC=ZBA*XDA-XBA*ZDA f 

AKLXAC=-YBA*XDA f 

DARG=AILXAC*AILXAC+AJLXAC*AJLXAC+AKLXAC*AKLXAC f 

ALXAC=DSQRT(DARG) 

C ALXAC= SQRTlAILXAC *A I LXAC+AJLXAC*AJLXAC+AKLXAC*AKLXAC ) M 

C*** the PRECEEDING YIELD THE UNIT NORMALS TO THE FLAT PLATE SEGMENTS fc 

C*** AREA DETERMINATION FOR TRAPEZOIDAL SEGMENT OF TWISTED FLAT PLATE $ 

DELX=CC( LI/NUMM1 

AREA=DELX*ALXAC/AC J’ 

C*** UNSTEADY PRESSURE FORCE ;pj 

FRCE~SUMP*AREA/ ALX AC II 

FXUS ( ItL ) = FRCE*AILXAC |I 

FYUS! I , L ) =FRCE*AJL X AC \ 

FZUS! I ,L)=FRCE*AKLXAC j 

c ********* RESULTANT VELOCITY OF BLADE CONTROL POINT RELATIVE TO ************ J 
C**************** BLADE-FIXED COORDINATE SYSTEM ********************* |j 

VXPOW=RPS s * c YCl L)-( WYC ( I ,L)«COSTH-WXC( I ,L)*SINTH) | 

VYPOW=-RPS*XC( I, L)-(WXCt I*L)#COSTH+W-YCt I ,L)*SINTH) 

VZPGW=-WZC( I, L) S _] 

Q********* ;*** ****************** *************** j 

PUS! I » L ) =F XUS ( I ,L)*VXPOW+FYUS( I , L ) *VYPOW+FZUS l I ,L)*VZPOW j 

P! I,L) = PQS ( I, L)+PUSCI*L) ? 

FXBV! I,L)=FXQS( I ,L )+FXUS(I,L) •! 

FYBV! I »L)=FYQS( I ,L ) +FYUS! I,L ) ‘ 

34 FZPV! I ,L )=FZOS( I ,L ) +FZUS! IiL) 

33 CONTINUE 
PRINT 1001 
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CP=G. 

C 400 

CT 1 = 0 . .. 

CP 1=0. " ,?■ 

PCWER=0. 

THRUST = 0 . 

DO 90 L = 1 • NO PAN 
THRST I L ) =0 . 

PG!;J!UL)=0. 

DRAG (*1)=0. 

TCRQ ( L ) =0 . 

C** SPANWISE DISTRIBUTION 
DO 91 1 = 1* NUMM1 
THRSTI L ) = THRST(L )+F2BV( l* L ) 

POWRIL)=PQWRlL )+P( I * Lrl 
DRAG(L)=DRAG(L) + FXBV( 1*1) 

91 TORQ(L)=DRAGIL)*YC (L) 

OLCT I P=Bl* THRSTI LI/COEF 

C DLCPIP=(BL*P0WRIL)/(C0EF«VTIP) )*PI**4/4. 

CT I=CT I +DL CT I P 
C CP I = CP I+DLCP I P 

DLCTIP=DLC TIP/I C Y ( L+ 1 > -Y ( L ) ) /R ) 

C CJLCPIP = DLCP IP/I (Y(L + 1)-Y(L) )/R) 

C PHI=ATAN(DRAG(L)/THRST(L) ) 

92 DARG=DRAG(L)/THRST I L ) 

PH I = DAT AN I OARG ) 

PRINT 1050* PHI 
DARG=PH I 

C 1 S1NFI=SIN(PHI ) 

C COSF I=COS (PHI) 

SINFI=DSIN(DARG) 

COSF1=DCOS(DARG) 

ALFA= I BETAC I L ) -PHI )*180»/PI 
IFITHCIL) .GE..08) A0=-0 ,0352*THC I L ) +0. 1 109 
IFITHCI L ) .GE. .21) A0=-0 . 1525*THC I L ) +0 . 138 1 5 
CL=AG*ALFA 

CDMIN=0.01563*THC( U+0.004 
C D=CDM I N 

VEVT=YCIL)*CQSFI/R 

SIGMA=BL*CC(L>/(PI*R) 

CED DELCTH=VEVT*VEVT*S I GMA* ( CL*COSF I -CD*S I NF I ) /2. 

CEO DELCPH=VEVT*VEVT*S I GMA* ( CL*S I NF I +CD*COSF I )*YC(L)/R/2. 

CED DELCTP=DELCTH*P I**3/4. 

CED DELCPP=DELCPH*PI**4/4. 

PRINT 1009, DLCTIP,DLCP I P » AL FA * DELCTP , DELC P P 
CED CT=CT+DELCTP*IY(L+ 1 ) - Y ( L ) )/R 
CED90 CP=CP+DELCPP*(Y(L+ 1 )-Y(L))/R 
90 CONTINUE 
PRINT 1002 

PRINT 1009, CTI , CP * CT ,CPI 
CED IF! ITIME.EQ.MXTIME ) GO TO 3004 
PRINT 1004 

C******************* JGAM FUR NEXT TIME STEP ***********#£******* 


CANADAIR 

CANADAIR 

CANADAIR 


nnoo oooo 


DG 41 I I = l f MATRX 1 
41 T G A M { II ) =G A KM A C 1 I ) 

C***************** SHED VORTICES ADDED ******************* 

L = Q n 

DO 42 I I =NUM f MATRX 1 1 NUM 
L = L + 1 

GAMS ( L » IT IMb+1) -GAMMA ( I I ) 

RASIL, ITIME+l > = GAMS(Li I T IME+1 ) / ( 2 . 0*P I * SUMARR ) 

4 2 RAS ( L i l T I ME+ l ) = ABS ( RAS ( L » ITIME + l) ) 

D0104 L= 1 » NOPAN 
104 PRINT1Q03»GAMS(L* 1 T IME+1 ) 

PRINT 1004 

$*$*********$* CONSERVATION OF ANGULAR M0MENTUM» SHED ***************** 
I F{ L INWA.NE. 1) GG TO 1195 
DO 43 L s 1 » NUPAN 

GAMS ( L t IT IME + 1 )=GAMS(L, ITIME+l )*SQRT ( ( X t NUM , 1+ 1 ) -X { NUM , L ) )**2+(Y<L 
1 + 1 )-Y(L) }**2 + ( L ( VUM » L+ 1 ) -7. ( NUM t L ) 1**2) 

43 CONTINUE 
1 19*5 CONTINUE 

DQ101 L= 1 « NOPAN 

101 PRINT 1003»GAMS(L» ITIME+'l) 

PR I NT 1004 

******************* STRENGTHS OF TRAILERS ADDED ****************** 

SUM 1 = 0. k 

DO 4 5 L= 1 1 NPANP 1 
M=(L-l)*NUM+l 
N=M+NUM-2 
SUM2=0. 

DO 44 I I = M » N 

IF(L.EQ.NOPAN+l)GO TO 44 
SUM2 = SUM2+GAMMA( II ) 

44 CONTINUE 

GAMT( L, ITIME+l )=SUM1-SUM2 

RAT (L, ITIME+l )=GAMT( L, I T IME+ 1 ) f ( 2 . 0*P I *SUMARR ) 

RAT ( L t IT IME+ l ) =ABS ( RAT ( L i ITIME+l l ) 

45 SUM1=SUM2 ! 

PRINT 7777 
DO 102 L= 1 » NPANP 1 

102 PRINT 8888,RAS(L, ITIME+l), RAT f L » I T IME+1 ) 

CE IF( ITIME.EQ.MXTIME )G0 TO 3004 

I F ( IT IMEA EQ.MXT I. ME ) GO TO 14 
C 102 PRINT1003,GAMT(L,I TIME+1) 

C PRINT10Q4 

C********** WAKE COORDINATE POSITION WRT PROPELLER DISC PLANE 
IT IMP1= IT IME+L 
DC 50 ITT= I » ITIMPI 
I T= I.T I ME-I TT +2 
DD SO L= 1 » NPANP 1 

C**** TRAILING EDGE SHED FILAMENT POSITIONS TRANSFORMED TO PROPELLER 
£*************«********* **C0ORDINATES 

XW ( L » 1 ) =X < NUM » L ) *C OSTH-Y ( L ) *51 NTH 

Y W ( L » 1 ) = Y ( L ) *C 0 S TH + X ( N UM » L ) * S I N T H — 

2WIL,1)=ZINUM,L ) 
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XD=XW { L » IT) 

YD=YW(L, IT) 

ZD S 7 W ( L » IT) 

£*************#**^^00 IT Y DUE TP BOUND VORTICITY ******************* 
VXB=0. 

VYB=0. 

VZB = 0 . 

I F( LINWA.EO. I ) GO TO 79 

COT = C.OSTH 

S I T = S I NTH 

H=HH 

I WAKE = 0 

CALL INVEL 

C I F I KMAX . EQ.NUM ) GO TO 1151 

C I F ( KMAX . EQ.NPANP 1 ) GO TO 1151 

C I F { ITT.GT. I) GO TO 1151 

C IF(L.GT.i) GO TO 1 151 

C PRINT 7777 

C DO 1137 K3= WNPANP 1 

C 1 137 PRINT 8888, ANIK3), RA(K3) 

C 1 137 PRINT 8888, RASIK3 t ITIME+1) ,RAT( K3 1 ITIME+1) 

1151 VXB=VXP 
VYB=VYP 
V ZB = VZ P 

79 CONTINUE 

c ************#veLOCITY AT WAKE POINTS DUE TO INTERACTION ***************** 
C**«=* WAKE COORDINATES IN PROPELLER AXIS SYSTEM 
VXW=0. 

V YW=0 . 

VZW=0. 

IFILINWA.EQ.1 ) GO TO 688 
I F ( ITIME.EQ.O) GO TO 688 
COT=l. 

SIT=0. 

H=HH 
I WAKE= 1 
CALL INVEL 

C I FIKMAX.EQ.NUM) GO TO 1152 

C IFIKMAX* EQ.NPANP 1 ) GO TO 1152 

C I F ( ITT.GT.l) GO TO 1152 

C IF(L-GT.l) GO TO 1152 

C PRINT 7777 

C DO 1135 K4= 1 , NP ANP 1 

C 1 1 35 PRINT 8888,AN(K4), RA(K4) 

C 1 135 PRINT 8868, RAS ( K4 , ITIME+1), RATIK4, ITIME+1) 

1152 VXW=VXP 
VYW=VYP 
VZW=VZP 

C******* LOCAL SELF-INDUCED VELOCITY ************* 

I TR A I L=0 

C 5GO 

X2=XQ 
Y2 = YD 
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Z2=ZD 

6B8 V1S=0. 

VJS=0. 

VKS=Q. 

I F( LINWA.EQ. L ) GO TO 68 
IF(ITIME.EQ.O) GO TO 66 
558 1 F ( ITRAID550, 550, 5 51 

550 IFIL.EQ.I.QR.L.EG. NPANP1) GO TO 552 
X 1 = XW ( L-l* IT) 

V 1- Y W ( L~ 1 , IT) 

Z 1= ZW ( L- 1 , IT) 

X3=XW ( L+l , IT) 

Y3=YW(L+1, IT) 

Z3=ZW l L+ 1 » I T ) 

DELS 2= SORT! { X3-X 2 ) **2+ ( Y3- Y2 ) **2+ I Z3- Z2 ) **2 ) 

DELS 1=SCRT ( (X2-X1) **2 + I Y2-Y 1 ) **2+ ( Z2-Z 1 ) **2 ) 

AK= (GAMS( L-l, ITT )/ DELS i+G AMS ( L » I TT ) /DELS2 ) /2. 

I F t LINWA.EQ. 1) GO TO 3000 
AK= l GAMS (L-l, ITTKGAMSa. ITT) ) /2 
CED RAS ( NPANP1 )=RAS ( NOPAN ) 

RASINPANPl, ITT )= RAS (NOP AN , ITT) 

3000 CRB=RAS(L, ITT) 

GC TG 553 

551 I F ( I T . EQ . IT IMP 1 ) GO TO 552 
£$***$**#$## £ND OF WAKE ^s##*####*#**#**** 

I F { IT.EQ.l ) GO TO 554 

£$$$#$£**£*#***#**#* TRAILING EDGE **#****#*#***$**#*tf * 
XL=XW{L, IT-l) 

Y1=YW(L, IT-1) 

Z1=ZW(L, IT-1) 

DELS 1= SORT ( IX2-X1) **2+( Y2-Y 1)**2+ (Z2~Z1)**2) 

AK1=GAMT(L, ITT)/DELS1 
IF(LINWA.EQ.l) GO TO 3001 
AK1=GAMT ( L, ITT ) 

C****** ITT+1 CORRESPONDS TO IT-1 ************** 

3001 GO TO 555 

554 Xl=X(NUMMltL)*COSTH-Y(L)*SINTH 
Y1=Y(L)*C0STH+X(NUMM1,L)*SINTH 
Z 1= Z ( NUMM1 1 L ) 

DELS 1=SGRT ( ( X2-X 1 ) **2+ ( Y2-Y 1 ) **2+ { ZZ-Z 1 ) **2 ) 

M=(L-1)*NUM+1 

NsM+NUM-2 

C*********** TRAILER STRENGTHS (+) FEEDING DOWNSTREAM ********** 
AK1 = 0 * 

DO 237 J=M,N 
IF(L-l) 238,239,238 

C ************ LEFT TIP TRAILER *************** 

239 AK1 = AK1-GAMMA( J ) 

GO TO 237 

238 IF(L-NPANPl) 241,240,241 

240 AK 1= AK 1+GAMMA ( J-NUM ) 

GO TO 237 

241 AK1=AK 1-GAMMA { J)+GAM'A( J-NUM) J 


onnnn 


A-17 


237 CONTINUE 
555 X3 = XW ( L , IT+1) 

Y3 = YW ( L , IT+1) 

Z3=ZW! L, IT+1 ) 

DELS2=SCRT( ( X3-X2 ) **2 + ( Y 3-Y2 ) **2 + ( Z3-Z2 ) **2 ) 
AK2=GAMT!L, ITT-l )/DELS2 
IFtLINWA.EQ. 1) GO TO 3002 
AK2=GAMT( L , ITT— l ) 

3002 CRB=RAT(L,ITT-i) 


GO TO 552 


AK={ AK1+AK2I/2. 

553 CONTINUE 

IFIDELSl.LT.CRBi OR. DELS2.LT. CRB) 

AK-AK*ALOG( l./CRB) / ( 4.*P I ) 

XX=( (X3-X?)/DELS2+!Xi-X2)/DELSi)/{ ( DELS 1+DELS2 ) /2 . ) 

Y Y= ( !Y3-Y2)/DELS2+(Y1-Y2)/DELS1)/! ( DELS 1 + DELS2 ) /2 . ) 

Z 2= ( ( Z3-Z2)/DELS2+ ( Z1-Z2 ) /DELS L ) / ( ( DE LS 1+DELS2 ) /2 . ) 

XXX=( (X3-X2 )*DELS1 /DELS2- { X1-X2 ) *DELS2/DELS 1 ) / ( DELS 1+DELS2 ) 

YYY= { ( Y3-Y2)*DELSl/D£LS2-( Y1-Y2 ) *DELS2/DELS1 ) / ( DELS 1+DELS2 ) 
ZZZ=( ( Z3-Z2 )*DELS1 /DELS2-! Z 1-Z2 ) *DELS2/DELS1 ) / (DELSI + DELS2 ) 

C IF! ITRAIL )2200, 2200, 2201 

C2200 PR I NT 200 5 

C GO TO 2202 

C2201 PRINT 2006 

C2202 CONTINUE 

2005 FORMAT ( * ','SHED SHED SHED SHED SHED SHED SGED SHED SHED SHED • 

2006 FORMAT ( 1 ', 'TRAIL TRAIL TRAIL TRAIL TRAIL TRAIL TRAIL TRAIL ') 

PRINT 2000, X3,Y3, Z3,DELS2 
PR I NT 200 I , X 1 , Y1 , Z 1 , DELS 1 
PRINT 2002, X2,Y2,Z2 
PRINT 2003, XD, YD, ZD 

AK.XX, YY 
, E15 . 8 , * 

, E15. 8 , * 

, E 15 . 8 , ' 

, E15.8, • 


2000 

2001 

2002 

2003 

2004 


FORMAT! 
FORMAT ( 
FORMAT! 
FORMAT! 
FORMAT! 


PR INT 2004, I TT, L 
•X3=' ,E15.8, » Y3 = 
•Xl=* ,£15.8, • Y 1= 
* X2= ' ,5 15.8, 


ZZ,XXX, YYY,ZZZ 


552 

556 


Y2 = 

, *XD=' ,E15.8, ' YD= 

,2! 14, 2X ) ,7! E14.5,2X) ) 
VIS=VIS+AK*( YYY*ZZ-ZZZ*YY) 
VJS=VJS+AK*(ZZZ*XX-XXX*ZZ) 
VKS=VKS+AK*!XXX*YY-YYY*XX ) 
CONTINUE 

IF! ITRAIL)556»556, 557 
ITRAIL**! 

PRINT 1004 


Z3= ' , E 1 5 • 8 , 
Zl= ' , E 1 5 . 8 , 
Z2= ' , E15. 8 ) 
ZD= ' , E 1 5 • 8 ) 


DELS2=* , E 1 5 . 8 ) 
DELS1=' , E 15 . 8 ) 


GO TO 558 
557 CONTINUE 

68 VI ( L , IT )=VXB+VXW+V+VIS 
V J ( L , IT )=VYB+VYW+V JS 
50 VK( L, IT )=VZ6+VZW+VKS 

INDUCED VELOCITIES AT WAKE POINTS WRT PROPELLER DISC PLANE 


PR 1 NT 1004 

IF(LINWA.EO.l) GO TO 74 
OC 72 IT= 1 , IT IMP 1 
CC 7 3 L= 1 , NPANP 1 
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C PRINT 1005, IT,XW(L, IT) ,YWIL, IT) ,ZW!L, IT) ,VI (tt IT) *VJ(L» IT) ,VK!L, IT) 

C 1 ,L 

D ( 1) -XW ( L , l T ) 

D( 2) =YW(L, IT) 

0(3) = Z W ( L , I T ) 

STATUS! 1 ) =0 

C CALL BLKWRIT15LTAPE1, 3, D, STATUS) 

73 CONTINUE 
72 CONTINUE 

STATUS! 1 )=1 

C CALL BLKWRIT! 5LTAPE1, 3, D, STATUS) 

PR INT1004 

7 A CONTINUE 

C ********CALCULATI0N OF WAKE COORDINATE POSITION ************* 

DO 69 L= 1 , NP ANP 1 
DO 69 I IT = I» IT IMP 1 
I T= I T IME - 1 1 T + 2 

XW!L, IT + 1 )=XW!L,IT )+VI(L, IT)*DELT 
YW ( L , IT + 1 ) =YW t L , IT ) +VJ ( L * I T ) *DELT 

69 ZW ( L , IT+1)=ZW(L, IT ) +VK ( L » I T ) *DELT 

C ************* CONSERVATION OF ANGULAR MOMENTUM , TRA I LERS *************** 
I F ( L INWA.NE. 1 ) GO TO 1196 
DO 70 L= 1 » NPANP 1 

GAMT ( L , ITIME+1 )-=GAMT(L, I T I M E+ 1 ) * SQR T < { XW < L , 2 ) -XW ( L , 1 ) ) **2+ ( YW ( L , 2 ) 
l-YW(Ltl) )**2+(ZW!L , 2 ) -ZW ( L , 1 ) )**2) 

70 CONTINUE 
C 600 

1196 CONTINUE 

IT IME=IT IMP1 
PRINT1014, ITIME 
GO TO 71 

C3004 PMIN=PM IN+100.0 
C E PRINT 5555 
CE PRINT 6666, PMIN 
CE I F ( PMIN. GT. 1000.0) GO TO 1192 
CE GO TO 3003 
14 CONTINUE 
C 1 1 9 l CONTINUE 
C CALLBLKREWD! 5LTAPE l ) 

1192 STOP 
END 

SUBROUTINE INVEL 
DOUBLE PRECISION DARG 

COMMON X ( 1 0 , 3 0 ) , Y ( 30 ) , Z ( 10 , 3 0 ) , XW < 3 0 , 1 00 ) , Y W ( 30 , 1 00 ) , ZW ! 30 , 1 

100) , GAMMA! 100 ) , GAM S ( 3 0 , 1 00 ) , GAM T ( 30 , 100 ) , AN ( 30 ) , RA ( 30 ) ,CMA(30) ,C! 
130 ) , BETA! 30), AS! 30 ) , U ( 30 , 31 ) , YYC ( 30 ) ,CC ( 30 ) ,BETAC!30) ,R AS! 30, 100) , 
1RAT130, 100 ) 

COMMON ~VXP , VYP , VZP , COT , S I T , I T I ME , I WAKE , BL ,7 3 L , NOP AN , NUN , XD , YD , ZD , 
1H,E,A0,R,KMAX,LINWA,V, SUMARR 
Pl=3. 1415927 
VXP=0. 

VYP=0. 

VZP=C. 


I F ( IWAKE) 1,1,2 
2 I TEST=Q 
GO TO 23 
1 ITEST=-1 
23 CONTINUE 

DO 7 I B= 1 , I 8L 
DARG=2.DO*PI*t IB-1 I/BL 
CSIN=DCCS(DARG) 

SS IN = DS IN( DARG ) 
SSIN=SIN{2.*PI*< IB-1 ) /BL ) 
CSIN=COS(2.*PI*[ IB-1J/BL) 
COSBL=COT-CS IN-SIT *SS IN 
SINBL=SIT*CSIN+COT*SSlN 
I F ( I T EST ) 4, 5 , 5 

4 JMAX=NOPAN 
KKAXaNUM 
KK=0 

GO TO 6 

5 JMAX=ITIME 

o DO 7 J= 1 » UMAX 
I F ( ITEST )8, 9, 9 
9 J1=J+1 

J2= JMAX-J+1 
KM AX=NOP AN 

I F { ITEST.GT.O ) KMA X=KMAX + 1 
8 DO 26 K= 1 , KMAX 
I F ( ITEST) 10c 11. 11 

10 JJ=KK*NUM+K 
GAM=GAMMA( jj ) 

CRA = H 

I F(K-NUM) 12, 13, 12 

12 K 1= 1 
K2 = 3 

GO TO 15 

13 K 1 = 2 
K2 = 2 

GO TO 15 

11 Kl=l 
K2=l 

15 DO 26 KKK=K1,K2 

I F ( ITEST ) 29,30,20 
29 GO TO ( 16,17,18) ,KKK 

16 XA=X (NUM, J)*COSEL-Y( J)*SINBl 
XB=X(K, J)*COSBL-Y{ J)*SINBL 
YA=Y ( J ) *CQSBL+X(NJK, J )*SINBL 
YB = Y( J)*COSBL+X(K, J )*SINBL 
ZA=Z ( NUM, J ) 

ZB=Z(K, J) 

GO TO 19 

17 XA=X(K, J1*C0SBU-Y( J)*SINBL 
XB=X(K, J+1)*C0SBL-Y( J+1)*SINBL 
YA=Y( J)*COSBL+X(K, J)*SINBL 

YB = Y ( J + 1)*C0SBL+X( K,J)*SINBL 


ZA=Z(K, J) 

ZB=Z(K, J+l) 

GO TO 19 

18 XA=X(K, J+i)*COSBL-Y( J+l )*SINBL 
XB = X (NUM, J+1)*CQ5BL-Y( J+i )*SINBL 
YA = Y( J + l )*CCSBL+X( K» J + L ) *S INBL 
YB=Y ( J+L )*COSBL + X( NUM , J + 1 ) #S I NBL 
ZA=Z ( K, J+L) 

ZB=Z(NUM, J+L) 

GO TO 19 
30 l_L = K 
LLL=K+1 
I I — J 1 
I 1 1 = J 1 

GAM^GAMS ( LLt J2 ) 

CED RASUL, J2)=GAMS(LL , J 2 ) ( { 2 . 0*P I *SUMARR ) 
CED RASUL, J2) = ABS(RAS (LL,J2) ) 

CRA=RAS ( LL, J2 ) 

GO TO 21 

20 LL=K 
LLL=K 
I I - J 

I I I = J 1 

GAM=GAMT(LL, J2) 

CED RAT( LL, J2)=GAMT(LL , J2)/( 2 . 0*P l *SUMARR ) 
CED RATUL, J2)=ABS(RAT (LL,J2) ) 

CRA=RAT ( LL » J2 ) 

21 XA=XKJLL, II )*COSBL-YW(LL, II)*SINBL. 

XB=X W (LLL , 1 1 1 )*C0S8L-YW(LLL, I I I )*SINBL 
YA = YW(LL» 1 1 )*COSBL+XW(LL, II )*SINBL 

Y B= YW ( LLL , I 1 1 )*COS BL+XW( LLL, I I I )*SINBL 
ZA = ZW( LL, 1 1 ) 

ZB= ZW (LLL » III ) 

19 XBA=XB-XA 
C 700 

YBA=YB-YA 
ZBA=ZB-ZA 
XDA=XD-XA 
YDA=YD-YA 
Z DA= ZD-Z A 
XOB=XD-XB 
YDB=YD”YB 
ZDB=ZD-ZB 

* ALS=XBA*XBA+YBA*YB A+ZBA*ZBA 
ACS=XDA*XDA+YDA*YD A+ZDA*ZDA 
BCS=XDB*XDB+YDB*YDB+ZDB*ZDB 
DARG= ALS 
AL=DSQRT ( DARG ) 

DARG= ACS 
AC=DSQRT ( DARG ) 

DARG=BCS 
BC=DSQRT( DARG) 

C AL= SORT ( ALS ) 


noononnnonoonnno 


C AC = SORT ( AC S ) 

C BC- SORT ( BCS ) 

A I LX AC=YB A*ZDA-ZEA *YDA 
AJLXAC = Z6A*X0A-ZDA*Xr*A 
AKLXAC=XBA*YDA-XDA*YBA 
1 1 AO I F ( IWAKE ) 31» 31, 34 
34 HH=H 

GO TO 32 

31 HH=E 

32 CONTINUE 

IFIAL.LT. CRA ) GO TO 26 

IFIAC.LT. CRA ) GO TO 26 

IFIBC.LT. CRA ) GO TO 26 

COSA= ( ACS+ALS-BCS) / ( AC*AL*2. ) 

T6MPA=ABS( l.-COSA*COSA) 

DARG=TEMPA 

HCrRE=AC*DSQRT( DARG > 

C HCORE=AC*SQRT(TEMPA) 

I F ( HCORE . LE.CRA ) GO TO 26 
CCSB= ( BCS+ALS-ACS) / ( BC*AL*2. ) 

IF(LINWA.NE.l) GO TO 24 
I F C IWAKE ) 24j 2A» 25 

23 AL= ACS 

24 VFN=GAM*(CQSA+COSB > / ( AL*ACS*TEMPA*P I *4 . ) 

VXP=VXP+VFN*AILXAC 

VYP= VYP+VFN^AJLXAC 
VZP=VZP+VFN*AKLXAC 

26 CONTINUE 

1148 I F ( ITEST)27,7,7 

27 KK= J 

7 CONTINUE 

I F ( ITEST ) 28f 22 t 28 
22 ITEST'-l 
GO TO 23 

28 RETURN 
END 

C /* PSUCC MXMLT 

jjt 

MXMLT 

SUBROUTINE TO MULTIPLY TWO MATRICES — SINGLE PRECISION 
A = VARIABLE NAME OF THE PREMULTIPLIER MATRIX 
B = VARIABLE NAME OF THE POSTMULT I P L I ER MATRIX 
C = VARIABLE NAME OF THE PRODUCT MATRIX 
M = NUMBER OF ROWS IN THE PREMULTIPLIER MATRIX 
N = NUMBER OF COLUMNS IN THE PREMULTIPLIER MATRIX 
K = NUMBER OF COLUMNS IN THE POSTMULTI PL I ER- MATRI X 
JA = NUMBER OF ROWS IN THE PREMULTIPLIER MATRIX AS DIMENSIONED 

JB = NUMBER OF ROWS IN THE POSTMULTIPLIER MATRIX AS DIMENSIONED 

JC = NUMBER OF ROWS IN THE PRODUCT MATRIX AS DIMENSIONED 
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TO FORTRAN IV/360 BY GANDER 
WATFDR COMPATIBILITY BY BOB GATSKI 

SUBROUTINE MXMLT (A , B, C» M, N, K, JA* JB, JC l 
DIMENSION A ( J A f M I, B(JB,K ), C ( JC , K ) 

DO 1 I - 1 » M 
DO L J= 1 » K 
SUM=0.0 
DO 2 L«1»N 

2 SUM = SUM + A (I ,L) *B( L, J ) 

1 C( I , J ) - SUM 
RETURN 
END 

SUBROUTINE MX I NV ( A , MD IM , N ) 

REAL A(MDIM,N) tBIGA, HOLD 
INTEGER L ( 100 ) , M ( 1 00 ) 

CONVERTED FROM SSP ROUTINE MINV BY R.S. BUTLER 

DO 80 K = l * N 
L(K)=K 
M ( K ) -K 
& I GA=A ( K f K ) 

DO 20 J”*K , N 
DO 20 I = K * N 

10 I F < A B S ( BIGAJ-ABSIA ( I , J ) ) ) 15,20,20 
15 BltaM I ,J ) 

L ( K j Pi 
M { K ) = J 
20 CONTINUE 
J=L < K ) 

IF(J-K) 35,35,25 
25 DO 30 1=1, N 
HQLD=-A { K , I ) 

A ( K , I ) = A ( J , I J 
30 A ( J , I ) =HOLD 
35 I =M ( K ) 

IF(I-K) 45,45,38 
38 DG 40 J = 1 » N 
HOLD=-A ( J , K ) 

A{ J,K )=A{ J, I) 

40 A ( J, I ) =HCLD 
45 DO 55 1=1, N 

IF(I-K) 50,55,50 
50 A ( I,K)=A( I , K ) / ( - B [ GA) 

55 CONTINUE 
DO 65 1=1, N 
HOLD= A ( I , K ) 

00 65 J = 1 » N 
IF(I-K) 60,65,60 
60 IFIJ-K) 62,65,62 
62 A( I , J )=FCLD*A(K, J) +A ( I , J ) 

65 CONTINUE 
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DQ 7 5 J=1,N 
IF(J-K) 70,75,70 
A(K,J )*A(K, JJ/BIGA 


00005250,' 


70 A(K,J)*A(K, JJ/BIGA 
BOO 

75 CONTINUE 

A ( K , K ) = 1 . 0/ B IGA 

SO CONTINUE 
K-N 

100 K=K-1 

I F ( K ) 150 , L50, 105 
105 I = L I K ) 

IF(I-K) 120,120,108 
108 DO 110 J= 1 , N 
HGLD=A( J,K) 

A ( J , K ) -“A ( J , IJ 
110 A C J , I )=HOLD 
120 J=M ( K ) 

IF(J-K) 100,100,125 
125 DO 130 1 = 1, N 
HCLD= A ( K , I ) 

AIK, I ) =-A ( J , I ) 

130 A( J, I ) =HOLO 
GO TO 100 
150 RETURN 
END 

SUBROUTINE GAUSSIN ) 

COMMON XXI 10,30) , Y ( 30 ) , Z ( 10 , 30 ) , XW ( 30 , 1 00 ) , Y W ( 30 , 1 00 ) , ZW ( 30 , ] 

100 ) , GAMMA! 100) , GAM S ( 30,100) ,GAMT(30, LOO ) , X(30) ,RA (30) , DMA (30) »Ci 
130) , BETA130) , AS( 30 ) , At 30, 31 ) ,YVC (30) ,CC( 30) ,BETAC(30) , RAS ( 30 , 100 ) , 
1RAT ( 30 , 100) 

COMMON VXP,VYP,VZP , COT, SIT, I TIME, IWAKE,BL, IB L , NOPAN , NUM , XD , YD , ZD , 
1H,E, AO,R,KMAX,LINWA,V, SUMARR 

SOLVE A SET OF N SIMULTANEOUS EQUATIONS WITH N UNKNOWNS BY USE 
OF GAUSSIAN ELIMINATION. ***NOTE*** ALWAYS INSURE THAT THE 

DIMENSION STATEMENT IS ALSO REGISTERED IN THE MAIN PROGRAM ***** 
TO SET UP THE PROGRAM THE NUMBER OF EQUATIONS IS N THE C65662 
ARE CALCULATED IN THE MAIN PROGRAM AND PLACED IN THE MATRIX A. 

THE PROGRAM SOLVES FOR THE UNKNOWNS X AND RETURNS TO THE MAI5 

PROGRAM 

NP=N+1 

NM=N-1 

DO 10 K= 1 , NM 

KP = K-U 

L=K 

DO 11 I =KP , N 

IF ( ABS ( A ( I,K) ).GT.ABS(A(L,K) ) ) L=I 
CONTINUE 

IF (A(L,L) ) 33 , 34 , 3 3 
34 PRINT, ’YOU CANT DO IT THIS WAY’ 

STOP 

33 IF (L.EQ.K) GO TO 71 
DO 70 J = K » NP 
TEMP= A ( K, J ) 


00005400! 


0000 5450. 


00005700 


00005750: 

00005800 

00005850 

00005900 

00005950 

00006000 


A(K'J)~A(LiJ) 

A I L , J J =TEMP 
?0 CONTINUE 
71 CONTINUE 

DO 10 I = K P » N 
B — A ( I f K ) / A ( K t K ) 

DO 10 J=KP»NP 
10 At I , J > = A t I * J )-B*A{ K, J ) 
X(N)=AIN,NP)/A(N,N ) 

DC 12 IN**1»NM 
I-N-IN 

xt n«AU »np) 

I P s I + 1 

DO L 3 J - I P » N 

13 Xt I ) = X ( I ) — Ai I, J)*X ( J) 
12 X ( I ) -X ( I ) / A ( 1,1) 

RETURN 

END 


// DATA . 

INPUT DU * 



20 

5 100 

1 


L * t ' 

0.0 

LO.O 

1.0 0.0 

0.025 

5.729578 

0.0 

0 . 333333330.0 

0.075 

5.729578 

0.0 

0 . 333333330.0 

0.125 

5 . 72957 B 

c.o 

0 . 333333330.0 

0 . 175 

5.729578 

0.0 

0 . 333333330.0 

0.225 

5.729570 

0.0 

0 . 333333330.0 

0.275 

5.729578 

0.0 

0 . 333333330.0 

0.325 

5.729578 

0.0 - 

0 . 333333330. Q 

0.375 

5.729578 

0.0 

0 . 333333330.0 

0 . A 25 

5.729578 

0.0 

0 . 333333330c 0 

0.475 

5.729578 

0.0 

0 . 333333330.0 

0.525 

5.729578 

0.0 

0 . 333333330.0 

0.575 

5.729578 

0.0 

0 . 333333330.0 

0.625 

5.729578 

0.0 

0 . 333333330.0 

0.675 

5.729578 

0.0 

0 . 333333330.0 

0.725 

5.729578 

0.0 

0 . 333333330.0 

0.775 

5.729578 

0.0 

0 . 333333330.0 

0.825 

5.729578 

0.0 

0 . 333333330.0 

0.875 

5.729578 

0.0 

0 . 333333330.0 

0.925 

5.729578 

0.0 

0 . 333333330.0 

0.975 

5.729578 

0.0 

0 . 333333330.0 

0.00 

5.729578 

0.0 

0 . 333333330.0 

0.05 

5 .7295 (8 

0.0 

0 . 333333330.0 

0.10 

5.729578 

0.0 

0 . 333333330.0 

0.15 

5.729578 

0.0 

0 . 333333330.0 

0.20 

5.729578 

0.0 

0 . 333333330.0 

0.25 

5.729578 

0.0 

0 . 333333330.0 

0.30 

5.729578 

0.0 

0 . 333333330.0 

0.35 

5.729578 

0.0 

0 . 333333330.0 

0.40 

5.729578 

0.0 

0 . 533333330.0 

0.45 

5.729578 

0.0 

0 . 333333330.0 

0.50 

5.729578 

0.0 

0 . 333333330.0 

0.55 

5.729578 

0.0 

0 . 333333330.0 

0.60 

5.729578 

0.0 

0 . 333333330.0 

0.65 

5.729578 

0.0 

0 . 333333330.0 

0.70 

5.729578 

0.0 

0 . 333333330.0 

0.75 

5.729578 

0.0 

0 . 333333330.0 

0.80 

5.729578 

0.0 

0 . 333333330.0 

0.85 

5.729578 

0.0 

0 . 333333330. 0 

0.90 

5.729578 

0.0 

0 . 333333330.0 

0.95 

5.729578 

0.0 

0 . 333333330.0 

1.00 

5.729578 

0.0 

0 . 333333330.0 




